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Abstract. A time-dependent search for neutrino 
flares from pre-defined directions in the whole sky 
is presented. The analysis uses a time clustering 
algorithm combined with an unbinned likelihood 
method. This algorithm provides a search for sig- 
nificant neutrino flares over time-scales that are 
not fixed a-priori and that are not triggered by 
multiwavelength observations. The event selection is 
optimized to maximize the discovery potential, taking 
into account different time-scales of source activity 
and background rates. The method is applied to a 
pre-defined list of bright and variable astrophysical 
sources using 22-string IceCube data. No significant 
excess is found. 
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I. Introduction 

Several astrophysical sources are known to have a 
variable photon flux at different wavelengths, showing 
flares that last between several minutes to several days. 
Hadronic models of Active Galactic Nuclei (AGNs) pre- 
dict 0][2] neutrino emission associated with these multi- 
wavelength (MWL) emissions. Time integrated analyses 
are less sensitive in this flaring scenario because they 
contain a higher background of atmospheric neutrinos 
and atmospheric muons. Therefore a time dependent 
analysis is more sensitive because it reduces the back- 
ground by searching smaller time scales around the flare. 
A direct approach that looks for this correlation using 
specific MWL observations is reported in Q. 

In order to make the flare search more general, and 
since MWL observations are scarce and not available for 
all sources, we take an approach not triggered by MWL 
observations. We apply a time-clustering algorithm 
(see H) to pre-defined source directions looking 
for the most significant accumulation in time (flare) 
of neutrino events over background, considering all 
possible combinations of event times. One disadvantage 
of this analysis is the increased number of trials, 
which reduces the significance. Nevertheless, for flares 
sufficiently shorter than the total observation period, 
the time clustering algorithm is more sensitive than a 
time integrated analysis. The predicted time scales are 
well below this threshold. 



II. Flare search algorithm 

The time clustering algorithm chooses the most 
promising flare time windows based on the times of the 



most signal-like events from the analyzed data. Each 
combination of these event times defines a search time 
window (Ati). For each Ati a significance parameter 
A, is calculated. The algorithm returns the best X max 
corresponding to the most significant cluster. The signif- 
icance can be obtained using two approaches: a binned 
method, as in the previous implementation |4|, and 
an improved unbinned maximum likelihood method [5 1 
which enhances the performance. 

The unbinned maximum likelihood method defines 
the significance parameter by: 



A = -2 log 
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where x s is the source location, h s and % are the best 
estimates of the number of signal events and source 
spectral index, respectively, which are found by max- 
imizing the likelihood, (£): 
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The background probability density function (pdf), 
Bi, calculated purely from data distributions, is given 
by: 

Bi = Pr ace (0i, 4H)Pt n ~' n (E i A)P? me (UA\ (3) 

where p^ pace describes the distribution of events in a 
given area (a zenith band of 8° is used for convenience). 
In a simple case this probability would be flat because of 
random distribution of background events. However, due 
to applied cuts, Earth absorption properties and detector 
geometry, this probability is dependent on zenith, 0i, and 
azimuth, <pi. The irregular azimuthal distribution caused 
by the detector geometry is shown in Fig. Q] For time 
integrated analyses covering one year the dependence 
on the azimuth is negligible because the exposure for 
all right ascension directions is integrated. However, an 
azimuth correction becomes important for time scales 
shorter than 1 day, reaching up to 40% difference, thus 
it should be included in time dependent analyses. p s P ace 
has value unity when integrated over solid angle inside 
the test region (i.e. zenith band). 

The energy probability p^ nergv is determined from 
the energy estimator distribution and depends on the 
zenith coordinate. In the southern sky an energy sensitive 
event selection is the most efficient way to reduce the at- 
mospheric muon background. This energy cut decreases 
with zenith angle, thus creating a zenith dependence 
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of the energy. Therefore a zenith dependent energy 
probability, shown in Fig. [2] is needed. Note that for 
the northern sky this correction is small. 
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Fig. 1 . Normalized azimuth distribution of the data sample reported 
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Fig. 2. Background energy pdf from data as a function of the energy 
estimator and zenith angle. The ultra high energy sample 1101 is used. 
The southern sky corresponds to cos(0) > 0. 

The time probability p* mie is described by fits to the 
event rates in the entire observed period as a function of 
time (ti). Two regions of the sky (South and North) are 
distinguished because they have different properties, thus 
the zenith, 8i, dependence of p^ zme . The northern sky 
sample consists mostly of atmospheric neutrinos which 
do not show a significant seasonal variation, therefore a 
constant fit is used. For the southern sky, a sinusoidal 
fit is used because it is dominated by a background of 
high energy atmospheric muons which have seasonal 
variation. These fits are shown in Fig. [3] and include the 
necessary correction for the uptimqj of the detector. It 
has been verified that the time modulations for different 
zenith bands within a half hemisphere are the same, thus 
allowing us to use all events inside the half hemisphere 
for the fit of the rates. 

'The uptime takes into account the inefficiency periods and data 
gaps after data quality selection. 




Fig. 3. Uptime corrected rates and their fits for the southern (left) 
and northern (right) skies. 



The signal pdf, Si, is given by: 

Si = P t space (\ Si - S s \,v i )Pr erav (E i) i ,'Y.)lt me , 

(4) 

where, the spatial probability, p^P ace j s a Gaussian func- 
tion of | xi — x s |, the space angular difference between 
the source location, x s , and each event's reconstructed 
direction, xi, and <7j, the angular error estimation of 
the reconstructed track. The estimator used for Oi is 
the size of the error ellipse around the maximum value 
of the reconstructed event track likelihood. The energy 
probability, p^ nergv ) constructed from signal simulation, 
is a function of the event energy estimation, Ei, the 
zenith coordinate, 9i, and the assumed energy spectral 
index of the source, 7 S (i? _7s ). A projection of p^ nergy 
for the whole sky is shown in Fig. [4] For a given 9i and 
j s the energy pdf is normalized to unity over Ei. For the 
energy a dedicated estimator of the number of photons 
per track length is used. The time probability, p? lme 5 
is constant since no flare time structure is assumed (i.e. 
taken to be flat in time). 
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Fig. 4. Projection for the whole sky of the energy component of the 
signal pdf as a function of the energy estimator and energy spectral 
index. The ultra high energy sample [ 1 1 is used. 

We use a binned method implementation of the time 
clustering algorithm as a crosscheck of our new un- 
binned analysis. In the case of the binned method, a 
circular angular search bin (2.5° radius) around the 
source direction is used. The times of the events that 
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(a) Southern sky at (dec=-7.6, ra=306.4) 
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(b) Northern sky at (dec=16.1, ra=343.5) 
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Fig. 5. Detection probability (3cr) for two source directions. The curves correspond to different time duration of the flares as function of the 
injected flux with a E~ 2 energy spectrum, using an unbinned time variable method (dashed), compared to a time integrated method (solid). 
The same mean number of events are injected into the time-windows (1, 7, 15, and 276 days) at each point on the x-axis, which is labeled 
with the equivalent flux corresponding to the full 276 day period. 



define the search time windows (At;) are given by all the 
events inside this angular bin. The significance parame- 
ter is obtained from Poisson statistics, given the number 
of expected background events inside the bin and the 
observed events in each cluster with multiplicity m. 
The expected number of background events is calculated 
by integrating, in the given time window, the fit to the 
rates, as described above. This calculation takes into 
account the zenith dependence of the background, in 
zenith bands with the size of the bin, the corresponding 
uptime factor and the azimuth correction. 

The best significance obtained for a cluster is cor- 
rected for trial factors by running several Monte Carlo 
background-only simulations. The simulation is done 
by creating distributions from data of zenith, azimuth, 
reconstruction error and energy estimator. The event 
characteristics are randomly taken from these distri- 
butions while considering the correlations between the 
different parameters. In order to study the performance 
of the algorithm, we calculate the neutrino flare detection 
probability as a function of the signal strength and 
duration of the flare by simulating signal events on top 
of background evento The properties of signal events 
are taken from a dedicated signal simulation and depend 
on the assumed energy spectral index. The Point Spread 
Function (PSF) is used to smear the events around 
the source location, thus simulating the effect of the 
direction reconstruction. For each simulation, a random 
time is chosen around which signal events are randomly 
injected inside the time window defined by the flare 
duration. The flare duration is investigated in the range 
from 1 day to 15 days, though the algorithm finds the 
best time window, which could be larger. We constrain 
the largest flare duration in the algorithm to be less than 

2 The integral of the Poisson distribution of the background events 
starts at (m-1) since the beginning and end of the time period are fixed 
from the data itself. 

3 The number of injected background and signal events is Poisson 
distributed. 



30 days, which is sensible from 7-ray observations. 

III. Source selection 

Since searching for all directions in the sky would de- 
crease the significance, we consider only a few promis- 
ing sources, thus reducing the number of trials. We 
select variable bright astrophysical sources in the whole 
sky. The selected blazars, including Flat Spectrum Radio 
Quasars (FSRQs) and Low-frequency peaked BL Lacs 
(LBLs), are taken from the confirmed Active Galactic 
Nuclei (AGN) in the third EGRET catalogue (3EG) 0. 
We also require that they are present in the current latest 
Fermi catalogue (0FGL) [7]. The criteria for selecting 
variable and bright source is based on the following 
parameters thresholds: 

. Variability index (3EG) > 1 

. Maximum 3EG flux (E > 100 MeV) > 40 [10~ 8 
ph cm~ 2 s -1 ] 

. Average 3EG flux (E > 100 MeV) > 15 [10~ 8 ph 

cm _2 s _1 ] 
• Inside visibility region of IceCube. 

The selected source list consists of 10 directions 
(Table |I| that are going to be tested with the time 
clustering algorithm. Models like [2] favor fluxes of 
higher energy neutrinos from FSRQ sources. Given the 
absorption of neutrinos at different energies in the Earth 
and the event cut strategy, southern sky FSRQs are more 
favored by these models because of their higher energy 
range of sensitivity. 

IV. Data Samples 

IceCubeO 22-string data from 2007-08 is used. It 
spans 310 days with an overall effective detector uptime 
of 88.9% (i.e. 276 days). The whole sky (declination 
range from -50° to 85°) is scanned. Different selection 
criteria are applied for the northern and southern skies. 
Previously obtained reconstructed datasets are used: the 
standard point source sample for the northern sky J9] 
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TABLE I 

Results for pre-defined variable astrophysical source candidates using the time clustering algorithm. 



Source 


Type 


dec [°] 


ra I°l 

ra l j 


p-vslue 




Fluence Limit (GeV/cm*^ ) 


GEV J0540-4359 


LBL 


-44.1 


84.7 


0.54 


7.08 


29.8 


GEV J 1626-2502 


FSRQ 


-25.5 


246.4 


0.41 


22.8 


22.8 


GEV J1832-2128 


FSRQ 


-21.1 


278.4 


0.64 


4.49 


12.0 


f^T^\T non 


FSRQ 


- /.O 


306.4 


1 


3.55 


3.7 


3C 279 


FSRQ 


-5.8 


194.1 


0.52 


0.19 


3.3 


3C 273 


FSRQ 


2.0 


187.3 


0.84 


1.97 


0.37 


CTA 102 


FSRQ 


11.7 


338.1 


0.27 


3.7 


1.42 


GEV J0530+1340 


FSRQ 


13.5 


82.7 


0.82 


3.48 


0.47 


3C 454.3 


FSRQ 


16.1 


343.5 


0.047 


0.5 


2.22 


GEV J0237+1648 


LBL 


16.6 


39.7 


0.59 


1.99 


1.08 



Note.- The flare duration of the best cluster is given by At. The fluence upper limit is calculated by integrating d^/dE X E over the 90% 

energy range and At, assuming a neutrino energy spectrum of E . 



(5114 events, declination from -5° to 85°, 1.4° sky- 
averaged median angular resolution) and the dedicated 
ultra high energy sample for the southern sky flOl (1877 
events in the whole sky, declination from -50° to 85°, 
1.3° sky-averaged median angular resolution). The first 
sample is optimized, within an unbinned method, for 
the optimal sensitivity to both hard and soft spectrum 
sources. The second sample was optimized for a binned 
method at ultra high energies. Therefore it should be 
noted that the binned method results are much better in 
the southern sky than in the northern sky. Nevertheless, 
the unbinned method, for an E~ 2 energy spectrum still 
performs better in the southern sky. 

The energy containment in these two regions is dif- 
ferent, with ranges from TeV to PeV and from PeV 
to EeV, in the northern and southern sky respectively. 
Event tracks are obtained with a multi-photoelectron □ 
(MPE) ifTTIl reconstruction which improves the angular 
resolution for high energies. 

V. Results 

The probability of a 3a flare detection using this time 
variable analysis (time clustering algorithm) for a given 
number of Poisson mean injected signal events with a 
E~ 2 energy spectrum is shown for two sources, at the 
southern and northern skies, for different time scales 
in Fig. |5] For comparison purposes, time-integrated 
detection probabilities integrated over the whole 22- 
string IceCube data period (276 days) are also given. 
For shorter flare durations the detection probability in- 
creases and is well above a time integrated search. In 
the northern sky, the same simulated signal of 5 mean 
injected events in a 7-day window was on average 1.3 
times more likely to be detected at 3cr with the unbinned 
time variable search than with the time integrated search, 
and in the southern sky, on average about twice as likely 
with the time variable search. 

The time clustering algorithm was applied to the 
selected sources candidates. No significant excess above 
the atmospheric background is found, therefore upper 

4 The MPE reconstruction takes the arrival time distribution of the 
first of N photons using the cumulative distribution of the single photon 
pdf. 



limits on the neutrino fluence are calculated. The results 
are presented in TableQ] The highest fluctuation observed 
corresponds to 3C 454.3 with a p-value of 4.7% (not 
including the trial factors due to looking at several 
sources). 

VI. Summary 

We have presented the sensitivity of the time cluster- 
ing algorithm using an unbinned maximum likelihood 
method. This is an improvement over the previous 
performances using a binned method and time integrated 
analyses. The search window for variable sources has 
been extended to the southern sky. IceCube 22-string 
data was analyzed using this method looking for neu- 
trinos flares with no a priori assumption on the time 
structure of the signal. Since no deviation from the 
background-only hypothesis was found, upper limits on 
the neutrino fluence from these sources were derived. 
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